%dXt = sigma(t)^2dt 
%sigma(t)=OU
clear all
% index J for X
betaK=15; IIpart=sparse(index2(betaK));
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
s=sigmahandle; a=1;b=1;m=1;sig0=1;
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=1+2*betaK+betaK*(betaK-1)/2; X0=zeros(row,1); 
% solver of X
disp('solving X')
tic
[time,X] = ode45(@randucousqdt,tspan,X0,options,T,a,b,m,sig0,betaK,s,IIpart);

%varitest(time,X.*X,T,'uncorrelated','oudt',num)

hold on
%%%the WCE price
price=zeros(51,1);
Npath=100000;
eta=randn(Npath,betaK);
etasquare=1/sqrt(2)*(eta.^2-1);
etaeta=zeros(Npath, betaK*(betaK-1)/2);
counter=1;
for ii=1:betaK-1
    for jj=ii+1:betaK
        etaeta(:,counter)=eta(:,ii).*eta(:,jj);
        counter=counter+1;
    end
end

for iter=1:51
asset=zeros(Npath,1);
for ii=2:betaK+1
    asset=asset+X(iter,ii)*eta(:,ii-1);
end
for ii=betaK+2:2*betaK+1
    asset=asset+X(iter,ii)*etasquare(:,ii-betaK-1);
end
for ii=2*betaK+2:2*betaK+1+betaK*(betaK-1)/2
    asset=asset+X(iter,ii)*etaeta(:,ii-2*betaK-1);
end
asset=asset+X(iter,1);

price(iter)=sum(abs(asset))/Npath;
end
toc
plot(time, price);

%%% the MC price
tic
intsigmasq=zeros(51,1);
for pn=1:Npath
    samplepath=getapath(5001,'ou');
    sigmasq=samplepath.^2; 
    for jj=1:51
        intsigmasq(jj)=intsigmasq(jj)+midptrule(sigmasq(1:(jj-1)*100+1),1/10000);
    end
end
toc
intsigmasq=intsigmasq/Npath;
plot(time,intsigmasq,'r')

